library(sva) # Surrogate Variable Adjust
library(tidyverse)
library(bladderbatch)
data(bladderdata)
#library(pamr)

pheno <- pData(bladderEset)

pheno |> glimpse()

edata <- exprs(bladderEset)

edata |> glimpse()

# full model include all variables
full.mod <- model.matrix(~ cancer, data = pheno)

# null model include only variable not of interest
# in this case we only have var of interest, set 1
mod0 <- model.matrix(~ 1, data = pheno)

# estimate number of latent factors ------
n.sv <- edata |>
  num.sv(full.mod, method = 'leek')

n.sv

# estimate value of latent factors
svobj <- edata |>
  sva(full.mod, mod0, n.sv = n.sv)

# when number of features > 100000
# filter for 2000 most variable features to speed up
n.sv <- edata |>
  num.sv(full.mod, method = 'leek', vfilter = 2000)

# estimate value of latent factors
svobj <- edata |>
  sva(full.mod, mod0, n.sv = n.sv, vfilter = 2000)

# perform F test unadjusted for latent factors----
qvalue <- edata |>
  f.pvalue(full.mod, mod0) |>
  p.adjust(method = 'BH')

table(qvalue < .05)

# perform latent-factors-adjusted F test ----
## include latent factors in full & null model
suro.var <- svobj$sv |> as_tibble()

full.modsv <- cbind(full.mod, suro.var)
mod0sv <- cbind(mod0, suro.var)

full.modsv |> head()
mod0sv |> head()

# f.pvalue need sv col not named?
qvalue.sv <- edata |>
  f.pvalue(full.modsv, mod0sv) |>
  p.adjust(method = 'BH')

# adjusted SDEG decrease from 15193 to 14813
table(qvalue.sv < .05)

# integrate with limma ---------
colnames(full.modsv) <- colnames(full.modsv) |>
  make.names()

limma.fit <- edata |>
  lmFit(full.modsv)

contrast.mat <- makeContrasts(contrasts = 'cancerCancer - cancerNormal',levels = full.modsv)

eb <- limma.fit |>
  contrasts.fit(contrast.mat) |>
  eBayes()

eb |>
  topTable(adjust.method = 'BH',number = Inf) |>
  filter(adj.P.Val < .05) |>
  nrow()

# ComBat for known batches ------
# supply batch vector, null model (-batch vec)
# adjusted edata object is returned
combat_edata <- edata |>
  ComBat(batch=pheno$batch,
         mod=mod0, par.prior=TRUE, 
         prior.plots=T)
# use param ebayes by default (par.prior=T)
# check prior plot to determine whether non-param (par.prior=F) setting need to use

# use mean.only=T if mild batch effect or variance difference among batches is expected
# use ref.batch if one batch is much bigger & better in quality

# ComBat for RNA-seq -------------
# use negative binominal regression
rna_mat <- rnbinom(400, size = 10, prob = .1) |>
  matrix(nrow = 50, ncol = 8)

batch <- c(rep(1,4), rep(2,4))

adjusted <- rna_mat |>
  ComBat_seq(batch = batch)

# specify group of interest
group <- rep(c(0,1), 4)

adjusted_group <- rna_mat |>
  ComBat_seq(batch = batch, group = group)

# can define multiple covariates
group2 <- c(0,0,1,1,0,0,1,1)

covar_mat <- cbind(group, group2)

adjusted_2group <- rna_mat |>
  ComBat_seq(batch = batch, covar_mod = covar_mat)

# remove known batch in linear model -----------
full.mod <- model.matrix(~ as.factor(cancer) + as.factor(batch), data = pheno)

null.mod <- model.matrix(~ as.factor(batch), data = pheno)

head(null.mod)

qvalue.lm <- edata |>
  f.pvalue(full.mod, null.mod) |>
  p.adjust(method = 'BH')

# fsva remove batch effects for prediction ----------
# need pamr lib
# freeze sv in train data to use in test data
set.seed(12345)

# split bladder into train data & test data
train.indicator <- sample(1:57, size = 30, replace = FALSE)
test.indicator <- (1:57)[-train.indicator]

train.data <- edata[train.indicator]
test.data <- edata[test.indicator]

train.pheno <- pheno[train.indicator,]
test.pheno <- pheno[test.indicator,]

mydata <- list(x = train.data, y = train.pheno$cancer)
mytrain <- pamr.train(mydata)

mytrain |>
  pamr.predict(test.data, threshold = 2) |>
  table(test.pheno$cancer)

# TODO

# svaseq for seq -------------
# use moderated log transform for counts data
bulkrna <- read_delim('mission/FPP/zww_sa_mice/data/himr_bulk.csv')

gene_abundant <- bulkrna |>
  filter(count > 5) |>
  count(gene_symbol) |>
  filter(n >= 2) |>
  pull(gene_symbol)

bulkrna <- bulkrna |>
  filter(gene_symbol %in% gene_abundant)

wide_bulk <- bulkrna |>
  pivot_wider(id_cols = -group, names_from = sample, values_from = count) |>
  column_to_rownames('gene_symbol')

group_bulk <- wide_bulk |>
  colnames() |>
  str_remove('.$') |>
  as.factor()

mod1 <- model.matrix(~group_bulk)

mod1
mod0 <- mod1[,1]

# n.sv can be auto determined
svaseq_bulk <- wide_bulk |>
  as.matrix() |>
  svaseq(mod1, mod0)

svaseq_bulk |> glimpse()

# pch plot character, like shape in ggplot2
svaseq_bulk$sv |> plot(pch = 19, col='blue')

# svaseq can be supervised with ERCC spike in ctrl
